**Title:** A multilevel adaptive reaction-splitting simulation method for stochastic reaction networks.

**Abstract: **In this work, we present a novel multilevel Monte Carlo method for kinetic simulation of stochastic reaction networks characterized by fast and slow reaction channels.

To produce efficient simulations, we automatically classify the reactions channels into the fast and slow classes. To this end, we first introduce the concept of the level of activity of a reaction channel, which depends on the current state of the system.

Then, we propose a low cost heuristic that allows us to adaptively split the set of reaction channels into two subsets characterized by either a high or low level of activity. Based on a time splitting technique, the increments associated with high activity channels are simulated using the tau-leap method while those associated with low activity channels are simulated using an exact method.

This path simulation technique, which we name mixed method, is amenable for coupled path generation and a corresponding multilevel Monte Carlo algorithm.

To estimate expected values of observables of the system at a prescribed final time, our method bounds the global computational error to be below a prescribed tolerance, TOL, within a given confidence level. This goal is achieved with a computational complexity of order O{TOL^{-2}}, the same as with a pathwise exact method, but with a smaller constant.

We also present a novel control variate technique based on the stochastic time change representation by Kurtz, which may dramatically reduce the variance of the coarsest level at a negligible computational cost.

Our numerical examples show substantial gains with respect to the standard Stochastic Simulation Algorithm (SSA).